home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MFSS.FOR < prev    next >
Encoding:
Text File  |  1985-11-29  |  6.6 KB  |  244 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MFSS
  5. C
  6. C        PURPOSE
  7. C           GIVEN A SYMMETRIC POSITIVE SEMI DEFINITE MATRIX , MFSS WILL
  8. C           (1) DETERMINE THE RANK AND LINEARLY INDEPENDENT ROWS AND
  9. C               COLUMNS
  10. C           (2) FACTOR A SYMMETRIC SUBMATRIX OF MAXIMAL RANK
  11. C           (3) EXPRESS NONBASIC ROWS IN TERMS OF BASIC ONES,
  12. C               EXPRESS NONBASIC COLUMNS IN TERMS OF BASIC ONES
  13. C               EXPRESS BASIC VARIABLES IN TERMS OF FREE ONES
  14. C           SUBROUTINE MFSS MAY BE USED AS A PREPARATORY STEP FOR THE
  15. C           CALCULATION OF THE LEAST SQUARES SOLUTION OF MINIMAL
  16. C           LENGTH OF A SYSTEM OF LINEAR EQUATIONS WITH SYMMETRIC
  17. C           POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX
  18. C
  19. C        USAGE
  20. C           CALL MFSS(A,N,EPS,IRANK,TRAC)
  21. C
  22. C        DESCRIPTION OF PARAMETERS
  23. C           A     - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC SEMI-
  24. C                   DEFINITE MATRIX STORED COLUMNWISE IN COMPRESSED FORM
  25. C                   ON RETURN A CONTAINS THE MATRIX T AND, IF IRANK IS
  26. C                   LESS THAN N, THE MATRICES U AND TU
  27. C           N     - DIMENSION OF GIVEN MATRIX A
  28. C           EPS   - TESTVALUE FOR ZERO AFFECTED BY ROUND-OFF NOISE
  29. C           IRANK - RESULTANT VARIABLE, CONTAINING THE RANK OF GIVEN
  30. C                   MATRIX A IF A IS SEMI-DEFINITE
  31. C                   IRANK = 0 MEANS A HAS NO POSITIVE DIAGONAL ELEMENT
  32. C                             AND/OR EPS IS NOT ABSOLUTELY LESS THAN ONE
  33. C                   IRANK =-1 MEANS DIMENSION N IS NOT POSITIVE
  34. C                   IRANK =-2 MEANS COMPLETE FAILURE, POSSIBLY DUE TO
  35. C                             INADEQUATE RELATIVE TOLERANCE EPS
  36. C           TRAC  - VECTOR OF DIMENSION N CONTAINING THE
  37. C                   SOURCE INDEX OF THE I-TH PIVOT ROW IN ITS I-TH
  38. C                   LOCATION, THIS MEANS THAT TRAC CONTAINS THE
  39. C                   PRODUCT REPRESENTATION OF THE PERMUTATION WHICH
  40. C                   IS APPLIED TO ROWS AND COLUMNS OF A IN TERMS OF
  41. C                   TRANSPOSITIONS
  42. C
  43. C        REMARKS
  44. C           EPS MUST BE ABSOLUTELY LESS THAN ONE. A SENSIBLE VALUE IS
  45. C           SOMEWHERE IN BETWEEN 10**(-4) AND 10**(-6)
  46. C           THE ABSOLUTE VALUE OF INPUT PARAMETER EPS IS USED AS
  47. C           RELATIVE TOLERANCE.
  48. C           IN ORDER TO PRESERVE SYMMETRY ONLY PIVOTING ALONG THE
  49. C           DIAGONAL IS BUILT IN.
  50. C           ALL PIVOTELEMENTS MUST BE GREATER THAN THE ABSOLUTE VALUE
  51. C           OF EPS TIMES ORIGINAL DIAGONAL ELEMENT
  52. C           OTHERWISE THEY ARE TREATED AS IF THEY WERE ZERO
  53. C           MATRIX A REMAINS UNCHANGED IF THE RESULTANT VALUE IRANK
  54. C           EQUALS ZERO
  55. C
  56. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  57. C           NONE
  58. C
  59. C        METHOD
  60. C           THE SQUARE ROOT METHOD WITH DIAGONAL PIVOTING IS USED FOR
  61. C           CALCULATION OF THE RIGHT HAND TRIANGULAR FACTOR.
  62. C           IN CASE OF AN ONLY SEMI-DEFINITE MATRIX THE SUBROUTINE
  63. C           RETURNS THE IRANK X IRANK UPPER TRIANGULAR FACTOR T OF A
  64. C           SUBMATRIX OF MAXIMAL RANK, THE IRANK X (N-IRANK) MATRIX U
  65. C           AND THE (N-IRANK) X (N-IRANK) UPPER TRIANGULAR TU SUCH
  66. C           THAT TRANSPOSE(TU)*TU=I+TRANSPOSE(U)*U
  67. C
  68. C     ..................................................................
  69. C
  70.       SUBROUTINE MFSS(A,N,EPS,IRANK,TRAC)
  71. C
  72. C
  73. C        DIMENSIONED DUMMY VARIABLES
  74.       DIMENSION A(1),TRAC(1)
  75.       DOUBLE PRECISION SUM
  76. C
  77. C        TEST OF SPECIFIED DIMENSION
  78.       IF(N)36,36,1
  79. C
  80. C        INITIALIZE TRIANGULAR FACTORIZATION
  81.     1 IRANK=0
  82.       ISUB=0
  83.       KPIV=0
  84.       J=0
  85.       PIV=0.
  86. C
  87. C        SEARCH FIRST PIVOT ELEMENT
  88.       DO 3 K=1,N
  89.       J=J+K
  90.       TRAC(K)=A(J)
  91.       IF(A(J)-PIV)3,3,2
  92.     2 PIV=A(J)
  93.       KSUB=J
  94.       KPIV=K
  95.     3 CONTINUE
  96. C
  97. C        START LOOP OVER ALL ROWS OF A
  98.       DO 32 I=1,N
  99.       ISUB=ISUB+I
  100.       IM1=I-1
  101.     4 KMI=KPIV-I
  102.       IF(KMI)35,9,5
  103. C
  104. C        PERFORM PARTIAL COLUMN INTERCHANGE
  105.     5 JI=KSUB-KMI
  106.       IDC=JI-ISUB
  107.       JJ=ISUB-IM1
  108.       DO 6 K=JJ,ISUB
  109.       KK=K+IDC
  110.       HOLD=A(K)
  111.       A(K)=A(KK)
  112.     6 A(KK)=HOLD
  113. C
  114. C        PERFORM PARTIAL ROW INTERCHANGE
  115.       KK=KSUB
  116.       DO 7 K=KPIV,N
  117.       II=KK-KMI
  118.       HOLD=A(KK)
  119.       A(KK)=A(II)
  120.       A(II)=HOLD
  121.     7 KK=KK+K
  122. C
  123. C        PERFORM REMAINING INTERCHANGE
  124.       JJ=KPIV-1
  125.       II=ISUB
  126.       DO 8 K=I,JJ
  127.       HOLD=A(II)
  128.       A(II)=A(JI)
  129.       A(JI)=HOLD
  130.       II=II+K
  131.     8 JI=JI+1
  132.     9 IF(IRANK)22,10,10
  133. C
  134. C        RECORD INTERCHANGE IN TRANSPOSITION VECTOR
  135.    10 TRAC(KPIV)=TRAC(I)
  136.       TRAC(I)=KPIV
  137. C
  138. C        MODIFY CURRENT PIVOT ROW
  139.       KK=IM1-IRANK
  140.       KMI=ISUB-KK
  141.       PIV=0.
  142.       IDC=IRANK+1
  143.       JI=ISUB-1
  144.       JK=KMI
  145.       JJ=ISUB-I
  146.       DO 19 K=I,N
  147.       SUM=0.D0
  148. C
  149. C        BUILD UP SCALAR PRODUCT IF NECESSARY
  150.       IF(KK)13,13,11
  151.    11 DO 12 J=KMI,JI
  152.       SUM=SUM-A(J)*A(JK)
  153.    12 JK=JK+1
  154.    13 JJ=JJ+K
  155.       IF(K-I)14,14,16
  156.    14 SUM=A(ISUB)+SUM
  157. C
  158. C        TEST RADICAND FOR LOSS OF SIGNIFICANCE
  159.       IF(SUM-ABS(A(ISUB)*EPS))20,20,15
  160.    15 A(ISUB)=DSQRT(SUM)
  161.       KPIV=I+1
  162.       GOTO 19
  163.    16 SUM=(A(JK)+SUM)/A(ISUB)
  164.       A(JK)=SUM
  165. C
  166. C        SEARCH FOR NEXT PIVOT ROW
  167.       IF(A(JJ))19,19,17
  168.    17 TRAC(K)=TRAC(K)-SUM*SUM
  169.       HOLD=TRAC(K)/A(JJ)
  170.       IF(PIV-HOLD)18,19,19
  171.    18 PIV=HOLD
  172.       KPIV=K
  173.       KSUB=JJ
  174.    19 JK=JJ+IDC
  175.       GOTO 32
  176. C
  177. C        CALCULATE MATRIX OF DEPENDENCIES U
  178.    20 IF(IRANK)21,21,37
  179.    21 IRANK=-1
  180.       GOTO 4
  181.    22 IRANK=IM1
  182.       II=ISUB-IRANK
  183.       JI=II
  184.       DO 26 K=1,IRANK
  185.       JI=JI-1
  186.       JK=ISUB-1
  187.       JJ=K-1
  188.       DO 26 J=I,N
  189.       IDC=IRANK
  190.       SUM=0.D0
  191.       KMI=JI
  192.       KK=JK
  193.       IF(JJ)25,25,23
  194.    23 DO 24 L=1,JJ
  195.       IDC=IDC-1
  196.       SUM=SUM-A(KMI)*A(KK)
  197.       KMI=KMI-IDC
  198.    24 KK=KK-1
  199.    25 A(KK)=(SUM+A(KK))/A(KMI)
  200.    26 JK=JK+J
  201. C
  202. C        CALCULATE I+TRANSPOSE(U)*U
  203.       JJ=ISUB-I
  204.       PIV=0.
  205.       KK=ISUB-1
  206.       DO 31 K=I,N
  207.       JJ=JJ+K
  208.       IDC=0
  209.       DO 28 J=K,N
  210.       SUM=0.D0
  211.       KMI=JJ+IDC
  212.       DO 27 L=II,KK
  213.       JK=L+IDC
  214.    27 SUM=SUM+A(L)*A(JK)
  215.       A(KMI)=SUM
  216.    28 IDC=IDC+J
  217.       A(JJ)=A(JJ)+1.D0
  218.       TRAC(K)=A(JJ)
  219. C
  220. C        SEARCH NEXT DIAGONAL ELEMENT
  221.       IF(PIV-A(JJ))29,30,30
  222.    29 KPIV=K
  223.       KSUB=JJ
  224.       PIV=A(JJ)
  225.    30 II=II+K
  226.       KK=KK+K
  227.    31 CONTINUE
  228.       GOTO 4
  229.    32 CONTINUE
  230.    33 IF(IRANK)35,34,35
  231.    34 IRANK=N
  232.    35 RETURN
  233. C
  234. C        ERROR RETURNS
  235. C
  236. C        RETURN IN CASE OF ILLEGAL DIMENSION
  237.    36 IRANK=-1
  238.       RETURN
  239. C
  240. C        INSTABLE FACTORIZATION OF I+TRANSPOSE(U)*U
  241.    37 IRANK=-2
  242.       RETURN
  243.       END
  244.